1 Background

We compare the inference results from TMB, aghq, adam, and tmbstan. Import these inference results as follows:

tmb <- readRDS("depends/tmb.rds")
aghq <- readRDS("depends/aghq.rds")
adam <- readRDS("depends/adam.rds")
tmbstan <- readRDS("depends/tmbstan.rds")

depends <- yaml::read_yaml("orderly.yml")$depends

Check that the parameters (latent field, hyperparameters, model outputs) sampled from each of the four methods are the same:

stopifnot(names(tmb$fit$sample) == names(aghq$quad$sample))
stopifnot(names(tmb$fit$sample) == names(adam$adam$sample))
stopifnot(names(tmb$fit$sample) == names(tmbstan$mcmc$sample))

1.1 Run details

For more information about the conditions under which these results were generated, see:

1.1.1 TMB

dependency_details <- function(i) {
  report_name <- names(depends[[i]])
  print(paste0("Inference results obtained from ", report_name, " with the query ", depends[[i]][[report_name]]$id))
  report_id <- orderly::orderly_search(query = depends[[i]][[report_name]]$id, report_name)
  print(paste0("Obtained report had ID ", report_id, " and was run with the following parameters:"))
  print(orderly::orderly_info(report_id, report_name)$parameters)
}

dependency_details(1)
## [1] "Inference results obtained from naomi-simple_fit with the query latest(parameter:tmb == TRUE)"
## [1] "Obtained report had ID 20230425-211024-cc2b7929 and was run with the following parameters:"
## $tmb
## [1] TRUE
## 
## $sample
## [1] TRUE
## 
## $aghq
## [1] FALSE
## 
## $k
## [1] 1
## 
## $grid_type
## [1] "product"
## 
## $s
## [1] 1
## 
## $tmbstan
## [1] FALSE
## 
## $hmc_laplace
## [1] FALSE
## 
## $niter
## [1] 1000
## 
## $nthin
## [1] 1
## 
## $adam
## [1] FALSE
## 
## $nsample
## [1] 1000
## 
## $area_level
## [1] 4

1.1.2 aghq

dependency_details(2)
## [1] "Inference results obtained from naomi-simple_fit with the query latest(parameter:aghq == TRUE && parameter:k == 3 && parameter:s == 8)"
## [1] "Obtained report had ID 20230424-192453-b044811f and was run with the following parameters:"
## $aghq
## [1] TRUE
## 
## $k
## [1] 3
## 
## $s
## [1] 8
## 
## $grid_type
## [1] "pca"
## 
## $sample
## [1] TRUE
## 
## $tmb
## [1] FALSE
## 
## $tmbstan
## [1] FALSE
## 
## $hmc_laplace
## [1] FALSE
## 
## $niter
## [1] 1000
## 
## $nthin
## [1] 1
## 
## $adam
## [1] FALSE
## 
## $nsample
## [1] 1000
## 
## $area_level
## [1] 4

1.1.3 adam

dependency_details(3)
## [1] "Inference results obtained from naomi-simple_fit with the query latest(parameter:adam == TRUE)"
## [1] "Obtained report had ID 20230425-190141-a3df2ca3 and was run with the following parameters:"
## $adam
## [1] TRUE
## 
## $sample
## [1] TRUE
## 
## $tmb
## [1] FALSE
## 
## $aghq
## [1] FALSE
## 
## $k
## [1] 1
## 
## $grid_type
## [1] "product"
## 
## $s
## [1] 1
## 
## $tmbstan
## [1] FALSE
## 
## $hmc_laplace
## [1] FALSE
## 
## $niter
## [1] 1000
## 
## $nthin
## [1] 1
## 
## $nsample
## [1] 1000
## 
## $area_level
## [1] 4

1.1.4 tmbstan

tmbstan_details <- dependency_details(4)
## [1] "Inference results obtained from naomi-simple_fit with the query latest(parameter:tmbstan == TRUE && parameter:niter > 50000)"
## [1] "Obtained report had ID 20230414-110613-89cd4244 and was run with the following parameters:"
## $tmbstan
## [1] TRUE
## 
## $niter
## [1] 1e+05
## 
## $nthin
## [1] 40
## 
## $tmb
## [1] FALSE
## 
## $aghq
## [1] FALSE
## 
## $k
## [1] 1
## 
## $grid_type
## [1] "product"
## 
## $s
## [1] 1
## 
## $hmc_laplace
## [1] FALSE
## 
## $adam
## [1] FALSE
## 
## $nsample
## [1] 5000
## 
## $area_level
## [1] 4
tmbstan_details
## $tmbstan
## [1] TRUE
## 
## $niter
## [1] 1e+05
## 
## $nthin
## [1] 40
## 
## $tmb
## [1] FALSE
## 
## $aghq
## [1] FALSE
## 
## $k
## [1] 1
## 
## $grid_type
## [1] "product"
## 
## $s
## [1] 1
## 
## $hmc_laplace
## [1] FALSE
## 
## $adam
## [1] FALSE
## 
## $nsample
## [1] 5000
## 
## $area_level
## [1] 4

1.2 Time taken

time_taken <- data.frame(
  "TMB" = tmb$time,
  "aghq" = aghq$time,
  "adam" = adam$time,
  "tmbstan" = tmbstan$time
)

write_csv(time_taken, "time_taken.csv")

time_taken
adam <- adam$adam

2 Histograms and ECDF difference plots

We create histograms and empirical cumulative distribution function (ECDF) difference plots of the samples from each method. All of the possible parameter names are as follows:

pars <- names(tmb$fit$sample)
pars
##  [1] "beta_alpha"                    "beta_anc_alpha"                "log_sigma_alpha_a"             "log_sigma_alpha_xa"            "u_alpha_xs"                    "log_sigma_rho_as"             
##  [7] "logit_phi_alpha_as"            "u_rho_xs"                      "u_rho_as"                      "anc_clients_t1_out"            "population_t1_out"             "plhiv_attend_t1_out"          
## [13] "anc_already_art_t1_out"        "anc_tested_pos_t1_out"         "anc_tested_neg_t1_out"         "ui_anc_rho_x"                  "log_sigma_ancrho_x"            "us_alpha_x"                   
## [19] "log_sigma_lambda_x"            "anc_rho_obs_t1_ll"             "hhs_prev_ll"                   "beta_anc_rho"                  "log_or_gamma"                  "u_alpha_a"                    
## [25] "log_sigma_rho_xa"              "log_sigma_or_gamma"            "u_rho_xa"                      "log_sigma_alpha_as"            "logit_phi_alpha_xs"            "logit_phi_rho_as"             
## [31] "us_rho_xs"                     "anc_art_new_t1_out"            "untreated_plhiv_attend_t1_out" "untreated_plhiv_num_t1_out"    "artattend_ij_t1_out"           "log_betaT"                    
## [37] "anc_known_pos_t1_out"          "artattend_t1_out"              "OmegaT_raw"                    "ui_anc_alpha_x"                "ui_lambda_x"                   "log_sigma_ancalpha_x"         
## [43] "u_rho_x"                       "u_alpha_x"                     "logit_phi_rho_x"               "us_rho_x"                      "beta_rho"                      "u_alpha_xa"                   
## [49] "logit_phi_alpha_a"             "u_rho_a"                       "logit_phi_rho_a"               "beta_lambda"                   "log_sigma_rho_a"               "u_alpha_as"                   
## [55] "us_alpha_xs"                   "log_sigma_alpha_xs"            "logit_phi_rho_xs"              "log_sigma_rho_xs"              "rho_t1_out"                    "anc_rho_t1_out"               
## [61] "artnum_t1_out"                 "anc_alpha_t1_out"              "alpha_t1_out"                  "infections_t1_out"             "anc_plhiv_t1_out"              "plhiv_t1_out"                 
## [67] "lambda_t1_out"                 "log_sigma_alpha_x"             "logit_phi_alpha_x"             "log_sigma_rho_x"               "hhs_artcov_ll"                 "artnum_t1_ll"

We will produce plots about the following subset of them. There is no particular reason to choose this subset rather than other, it’s quite arbitrary.

pars_eval <- pars %in% c("beta_rho", "beta_alpha", "beta_lambda", "beta_anc_rho", "beta_anc_alpha", "u_rho_x", "us_rho_x", "u_rho_xs")
names(pars_eval) <- pars

We will write the tmbstan \(\hat R\) and ESS for the particular parameter on each plot to help in deciding how seriously to take the MCMC results:

rhats <- bayesplot::rhat(tmbstan$mcmc$stanfit)
ess_ratio <- bayesplot::neff_ratio(tmbstan$mcmc$stanfit)
niter <- 4 * tmbstan_details$niter / tmbstan_details$nthin
ess <- ess_ratio * niter

The beta_rho and beta_anc_rho parameters are used in other reports, so create dataframes here to save as artefacts:

beta_rho <- histogram_and_ecdf("beta_rho", i = 1, return_df = TRUE)
saveRDS(beta_rho$df, "beta_rho.rds")

beta_anc_rho <- histogram_and_ecdf("beta_anc_rho", return_df = TRUE)
saveRDS(beta_anc_rho$df, "beta_anc_rho.rds")

2.1 beta_rho

histogram_and_ecdf_helper <- function(par) lapply(1:sum(names(tmb$fit$obj$env$par) == par), histogram_and_ecdf, par = par)
histogram_and_ecdf_helper("beta_rho")
## [[1]]

## 
## [[2]]

2.2 beta_alpha

histogram_and_ecdf_helper("beta_alpha")
## [[1]]

## 
## [[2]]

2.3 beta_lambda

histogram_and_ecdf_helper("beta_lambda")
## [[1]]

## 
## [[2]]

2.4 beta_anc_rho

histogram_and_ecdf("beta_anc_rho")

2.5 beta_anc_alpha

histogram_and_ecdf("beta_anc_alpha")

2.6 logit

lapply(pars[stringr::str_starts(pars, "logit")], histogram_and_ecdf)

2.7 log_sigma

lapply(pars[stringr::str_starts(pars, "log_sigma")], histogram_and_ecdf)

2.8 u_rho_x

histogram_and_ecdf_helper("u_rho_x")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

## 
## [[18]]

## 
## [[19]]

## 
## [[20]]

## 
## [[21]]

## 
## [[22]]

## 
## [[23]]

## 
## [[24]]

## 
## [[25]]

## 
## [[26]]

## 
## [[27]]

## 
## [[28]]

## 
## [[29]]

## 
## [[30]]

## 
## [[31]]

## 
## [[32]]

2.9 us_rho_x

histogram_and_ecdf_helper("us_rho_x")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

## 
## [[18]]

## 
## [[19]]

## 
## [[20]]

## 
## [[21]]

## 
## [[22]]

## 
## [[23]]

## 
## [[24]]

## 
## [[25]]

## 
## [[26]]

## 
## [[27]]

## 
## [[28]]

## 
## [[29]]

## 
## [[30]]

## 
## [[31]]

## 
## [[32]]

2.10 u_rho_xs

histogram_and_ecdf_helper("u_rho_xs")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

## 
## [[18]]

## 
## [[19]]

## 
## [[20]]

## 
## [[21]]

## 
## [[22]]

## 
## [[23]]

## 
## [[24]]

## 
## [[25]]

## 
## [[26]]

## 
## [[27]]

## 
## [[28]]

## 
## [[29]]

## 
## [[30]]

## 
## [[31]]

## 
## [[32]]

2.11 us_rho_xs

histogram_and_ecdf_helper("us_rho_xs")

2.12 u_rho_as

histogram_and_ecdf_helper("u_rho_as")

2.13 u_alpha_x

histogram_and_ecdf_helper("u_alpha_x")

2.14 us_alpha_x

histogram_and_ecdf_helper("us_alpha_x")

2.15 u_alpha_xs

histogram_and_ecdf_helper("u_alpha_xs")

2.16 us_alpha_xs

histogram_and_ecdf_helper("us_alpha_xs")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

## 
## [[18]]

## 
## [[19]]

## 
## [[20]]

## 
## [[21]]

## 
## [[22]]

## 
## [[23]]

## 
## [[24]]

## 
## [[25]]

## 
## [[26]]

## 
## [[27]]

## 
## [[28]]

## 
## [[29]]

## 
## [[30]]

## 
## [[31]]

## 
## [[32]]

2.17 u_alpha_a

histogram_and_ecdf_helper("u_alpha_a")

2.18 u_alpha_as

histogram_and_ecdf_helper("u_alpha_as")

2.19 u_alpha_xa

histogram_and_ecdf_helper("u_alpha_xa")

2.20 ui_lambda_x

histogram_and_ecdf_helper("ui_lambda_x")

2.21 ui_anc_rho_x

histogram_and_ecdf_helper("ui_anc_rho_x")

2.22 ui_anc_alpha_x

histogram_and_ecdf_helper("ui_anc_alpha_x")

2.23 log_or_gamma

histogram_and_ecdf_helper("log_or_gamma")

2.24 Outcome variables

These are all of the outcome variables which we might be interested in to monitor:

names(tmb$fit$sample)[!(names(tmb$fit$sample) %in% unique(names(tmb$fit$obj$env$par)))]
##  [1] "anc_clients_t1_out"            "population_t1_out"             "plhiv_attend_t1_out"           "anc_already_art_t1_out"        "anc_tested_pos_t1_out"         "anc_tested_neg_t1_out"        
##  [7] "anc_rho_obs_t1_ll"             "hhs_prev_ll"                   "u_rho_xa"                      "anc_art_new_t1_out"            "untreated_plhiv_attend_t1_out" "untreated_plhiv_num_t1_out"   
## [13] "artattend_ij_t1_out"           "anc_known_pos_t1_out"          "artattend_t1_out"              "rho_t1_out"                    "anc_rho_t1_out"                "artnum_t1_out"                
## [19] "anc_alpha_t1_out"              "alpha_t1_out"                  "infections_t1_out"             "anc_plhiv_t1_out"              "plhiv_t1_out"                  "lambda_t1_out"                
## [25] "hhs_artcov_ll"                 "artnum_t1_ll"

The most important ones (this has been checked with Jeff) are HIV prevalence (rho_t1_out), ART coverage (alpha_t1_out), and HIV incidence (lambda_t1_out).

  • HIV prevalence has 6417 variables (rows)
  • ART coverage has 6417 variables (rows)
  • HIV incidence has 6417 variables (rows)
dim(tmb$fit$sample$rho_t1_out)
## [1] 6417 1000
dim(tmb$fit$sample$alpha_t1_out)
## [1] 6417 1000
dim(tmb$fit$sample$lambda_t1_out)
## [1] 6417 1000

Within tmb$naomi_data there is a dataframe called mf_out (modelframe output) which contains the area, age, sex mapping these 6417 rows. As a sanity check, confirm that the number of areas times number of ages times number of sexes indeed equals the number of rows:

mf_out <- tmb$naomi_data$mf_out
(n_area <- length(unique(mf_out$area_id)))
## [1] 69
(n_age <- length(unique(mf_out$age_group)))
## [1] 31
(n_sex <- length(unique(mf_out$sex)))
## [1] 3
stopifnot(n_area * n_age * n_sex == dim(tmb$fit$sample$rho_t1_out)[1])

It makes sense only to assess inferential accuracy at the finest level so as to avoid double counting.

mf_out$area_id %>% unique()
##  [1] "MWI"           "MWI_1_1_demo"  "MWI_1_2_demo"  "MWI_1_3_demo"  "MWI_2_1_demo"  "MWI_2_2_demo"  "MWI_2_3_demo"  "MWI_2_4_demo"  "MWI_2_5_demo"  "MWI_3_1_demo"  "MWI_3_10_demo" "MWI_3_11_demo"
## [13] "MWI_3_12_demo" "MWI_3_13_demo" "MWI_3_14_demo" "MWI_3_15_demo" "MWI_3_16_demo" "MWI_3_17_demo" "MWI_3_18_demo" "MWI_3_19_demo" "MWI_3_2_demo"  "MWI_3_20_demo" "MWI_3_21_demo" "MWI_3_22_demo"
## [25] "MWI_3_23_demo" "MWI_3_24_demo" "MWI_3_25_demo" "MWI_3_26_demo" "MWI_3_27_demo" "MWI_3_28_demo" "MWI_3_3_demo"  "MWI_3_4_demo"  "MWI_3_5_demo"  "MWI_3_6_demo"  "MWI_3_7_demo"  "MWI_3_8_demo" 
## [37] "MWI_3_9_demo"  "MWI_4_1_demo"  "MWI_4_10_demo" "MWI_4_11_demo" "MWI_4_12_demo" "MWI_4_13_demo" "MWI_4_14_demo" "MWI_4_15_demo" "MWI_4_16_demo" "MWI_4_17_demo" "MWI_4_18_demo" "MWI_4_19_demo"
## [49] "MWI_4_2_demo"  "MWI_4_20_demo" "MWI_4_21_demo" "MWI_4_22_demo" "MWI_4_23_demo" "MWI_4_24_demo" "MWI_4_25_demo" "MWI_4_26_demo" "MWI_4_27_demo" "MWI_4_28_demo" "MWI_4_29_demo" "MWI_4_3_demo" 
## [61] "MWI_4_30_demo" "MWI_4_31_demo" "MWI_4_32_demo" "MWI_4_4_demo"  "MWI_4_5_demo"  "MWI_4_6_demo"  "MWI_4_7_demo"  "MWI_4_8_demo"  "MWI_4_9_demo"
mf_out$sex %>% unique()
## [1] "both"   "female" "male"
mf_out$age_group %>% unique()
##  [1] "Y000_004" "Y000_014" "Y000_064" "Y000_999" "Y005_009" "Y010_014" "Y010_019" "Y015_019" "Y015_024" "Y015_049" "Y015_064" "Y015_999" "Y020_024" "Y025_029" "Y025_034" "Y025_049" "Y030_034" "Y035_039"
## [19] "Y035_049" "Y040_044" "Y045_049" "Y050_054" "Y050_064" "Y050_999" "Y055_059" "Y060_064" "Y065_069" "Y065_999" "Y070_074" "Y075_079" "Y080_999"
mf_out_fine <- mf_out %>%
  tibble::rownames_to_column("id") %>%
  mutate(id = as.numeric(id)) %>%
  filter(
    area_id %in% paste0("MWI_4_", 1:32, "_demo"),
    sex %in% c("male", "female"),
    age_group %in%
      c(
        "Y000_004", "Y005_009", "Y010_014", "Y015_019", "Y020_024", "Y025_029",
        "Y025_034", "Y030_034", "Y035_039", "Y040_044", "Y045_049", "Y050_054",
        "Y055_059", "Y060_064", "Y065_069", "Y070_074", "Y075_079", "Y080_999"
      )
  )

An id column has been added to mf_out_fine to enable merging to the samples.

head(mf_out_fine)
nrow(mf_out_fine)
## [1] 1152
rho_t1_out_fine <- tmb$fit$sample$rho_t1_out[mf_out_fine$id, ]
dim(rho_t1_out_fine)
## [1] 1152 1000
rho_t1_out_fine[1:10, 1:10]
##              [,1]        [,2]        [,3]        [,4]        [,5]        [,6]        [,7]        [,8]        [,9]       [,10]
##  [1,] 0.003475839 0.003667526 0.003772017 0.003670830 0.003084421 0.003381271 0.003227835 0.003240394 0.003243281 0.003140848
##  [2,] 0.005383996 0.005680915 0.005842769 0.005686032 0.004777698 0.005237512 0.004999844 0.005019297 0.005023768 0.004865102
##  [3,] 0.006455132 0.006811122 0.007005176 0.006817257 0.005728212 0.006279505 0.005994553 0.006017877 0.006023238 0.005833005
##  [4,] 0.010537810 0.013036405 0.010803304 0.012115244 0.008440319 0.009413301 0.009728890 0.009145472 0.009130840 0.008276510
##  [5,] 0.027403103 0.026320700 0.028744722 0.027666183 0.020696399 0.025604750 0.023652769 0.022730970 0.024487682 0.026568189
##  [6,] 0.053177866 0.052120950 0.056024483 0.049817320 0.043130748 0.045152278 0.045047404 0.048507542 0.045410425 0.047294740
##  [7,] 0.062490740 0.066736508 0.069409553 0.063366683 0.054919044 0.057670696 0.058735515 0.057784846 0.056701005 0.055864714
##  [8,] 0.075608754 0.087323825 0.088263617 0.082452168 0.071523909 0.075304006 0.078016440 0.070852756 0.072604794 0.067936286
##  [9,] 0.106334220 0.112370623 0.111696858 0.113255673 0.099951306 0.105550418 0.094606976 0.097801372 0.107369936 0.103324468
## [10,] 0.114587139 0.110216218 0.113231040 0.119736201 0.101027883 0.117628993 0.099670952 0.107664312 0.100897850 0.098092128

Something weird happening with adam, go back and check this:

histogram_and_ecdf(par = "alpha_t1_out", i = 1)

3 KS plots

3.1 Individual parameters

ks_helper <- function(par, starts_with = FALSE, ...) {
  to_ks_df(par, starts_with = starts_with) %>% ks_plot(par, ...)
}

3.1.1 beta

ks_helper("beta", starts_with = TRUE)

3.1.2 logit

ks_helper("logit", starts_with = TRUE)
## Warning: Removed 116 rows containing missing values (`geom_tile()`).

3.1.3 log_sigma

ks_helper("log_sigma", starts_with = TRUE)
## Warning: Removed 116 rows containing missing values (`geom_tile()`).

3.1.4 u_rho_x

ks_helper("u_rho_x")

3.1.5 u_rho_xs

ks_helper("u_rho_xs")

3.1.6 us_rho_x

ks_helper("us_rho_x")

3.1.7 us_rho_xs

ks_helper("us_rho_xs")

3.1.8 u_rho_a

ks_helper("u_rho_a")

3.1.9 u_rho_as

ks_helper("u_rho_as")

3.1.10 u_alpha_x

ks_helper("u_alpha_x")

3.1.11 u_alpha_xs

ks_helper("u_alpha_xs")

3.1.12 us_alpha_x

ks_helper("us_alpha_x")

3.1.13 us_alpha_xs

ks_helper("us_alpha_xs")

3.1.14 u_alpha_a

ks_helper("u_alpha_a")

3.1.15 u_alpha_as

ks_helper("u_alpha_as")

3.1.16 u_alpha_xa

ks_helper("u_alpha_xa")

3.1.17 ui_anc_rho_x

ks_helper("ui_anc_rho_x")

3.1.18 ui_anc_alpha_x

ks_helper("ui_anc_alpha_x")

3.1.19 log_or_gamma

ks_helper("log_or_gamma")

3.1.20 Outcome variables

ks_df_out <- function(par) {
  to_ks_df(par = par, outputs = TRUE, id = mf_out_fine$id) %>%
  select(-par) %>%
  left_join(
    mf_out_fine %>%
      rename("full_id" = "id") %>%
      tibble::rowid_to_column("index") %>%
      mutate(index = as.numeric(index))
  )
}

ks_rho_t1_out <- ks_df_out(par = "rho_t1_out")
## Joining, by = "index"
ks_alpha_t1_out <- ks_df_out(par = "alpha_t1_out")
## Joining, by = "index"
ks_lambda_t1_out <- ks_df_out(par = "lambda_t1_out")
## Joining, by = "index"
ks_plot(ks_df = ks_rho_t1_out, par = "rho_t1_out", alpha = 0.2)

ks_plot(ks_df = ks_alpha_t1_out, par = "alpha_t1_out", alpha = 0.2)

ks_plot(ks_df = ks_lambda_t1_out, par = "lambda_t1_out", alpha = 0.2)
## Warning: Removed 116 rows containing missing values (`geom_tile()`).

3.2 Summary table

Values in dark blue are the lowest KS difference for that particular parameter:

options(dplyr.summarise.inform = FALSE)

ks_summary <- lapply(unique(names(tmb$fit$obj$env$par)), function(x) {
  to_ks_df(x) %>%
    group_by(method, index) %>%
    summarise(ks = mean(ks), par = x) %>%
    ungroup()
}) %>%
  bind_rows() %>%
  pivot_wider(names_from = "method", values_from = "ks") %>%
  rename(
    "Parameter" = "par",
    "KS(adam, tmbstan)" = "adam",
    "KS(aghq, tmbstan)" = "aghq",
    "KS(TMB, tmbstan)" = "TMB",
  )

r <- adam$quad$obj$env$random
x_names <- names(adam$quad$obj$env$par[r])
theta_names <- names(adam$quad$obj$env$par[-r])

dict <- data.frame(
  Parameter = c(unique(x_names), unique(theta_names)),
  Type = c(rep("Latent field", length(unique(x_names))), rep("Hyper", length(unique(theta_names))))
)

ks_summary <- ks_summary %>%
  left_join(dict, by = "Parameter")

ks_summary %>%
  gt::gt() %>%
  gt::fmt_number(
    columns = starts_with("KS"),
    decimals = 3
  ) %>%
  gt::tab_style(
    style = cell_fill(color = "#0072B2"),
    locations = cells_body(
      columns = `KS(adam, tmbstan)`,
      rows = `KS(adam, tmbstan)` < `KS(aghq, tmbstan)` & `KS(adam, tmbstan)` < `KS(TMB, tmbstan)`
    )
  ) %>%
    gt::tab_style(
    style = cell_fill(color = "#0072B2"),
    locations = cells_body(
      columns = `KS(aghq, tmbstan)`,
      rows = `KS(aghq, tmbstan)` < `KS(adam, tmbstan)` & `KS(aghq, tmbstan)` < `KS(TMB, tmbstan)`
    )
  ) %>%
    gt::tab_style(
    style = cell_fill(color = "#0072B2"),
    locations = cells_body(
      columns = `KS(TMB, tmbstan)`,
      rows = `KS(TMB, tmbstan)` < `KS(aghq, tmbstan)` & `KS(TMB, tmbstan)` < `KS(adam, tmbstan)`
    )
  )
index Parameter KS(adam, tmbstan) KS(aghq, tmbstan) KS(TMB, tmbstan) Type
1 beta_rho 0.073 0.080 0.076 Latent field
2 beta_rho 0.074 0.065 0.074 Latent field
1 beta_alpha 0.069 0.040 0.068 Latent field
2 beta_alpha 0.111 0.111 0.102 Latent field
1 beta_lambda 0.031 0.095 0.147 Latent field
2 beta_lambda 0.023 0.046 0.043 Latent field
1 beta_anc_rho 0.032 0.113 0.092 Latent field
1 beta_anc_alpha 0.084 0.077 0.046 Latent field
1 logit_phi_rho_x 0.562 0.562 0.139 Hyper
1 log_sigma_rho_x 0.642 0.642 0.157 Hyper
1 logit_phi_rho_xs 0.535 0.535 0.104 Hyper
1 log_sigma_rho_xs 0.669 0.669 0.277 Hyper
1 logit_phi_rho_a 0.508 0.508 0.065 Hyper
1 log_sigma_rho_a 0.566 0.566 0.105 Hyper
1 logit_phi_rho_as 0.545 0.545 0.108 Hyper
1 log_sigma_rho_as 0.570 0.570 0.152 Hyper
1 log_sigma_rho_xa 0.680 0.680 0.207 Hyper
1 u_rho_x 0.050 0.055 0.082 Latent field
2 u_rho_x 0.041 0.105 0.075 Latent field
3 u_rho_x 0.031 0.020 0.035 Latent field
4 u_rho_x 0.037 0.090 0.090 Latent field
5 u_rho_x 0.125 0.062 0.091 Latent field
6 u_rho_x 0.125 0.152 0.133 Latent field
7 u_rho_x 0.084 0.157 0.174 Latent field
8 u_rho_x 0.059 0.041 0.045 Latent field
9 u_rho_x 0.042 0.067 0.052 Latent field
10 u_rho_x 0.052 0.082 0.125 Latent field
11 u_rho_x 0.062 0.045 0.035 Latent field
12 u_rho_x 0.116 0.149 0.175 Latent field
13 u_rho_x 0.038 0.050 0.039 Latent field
14 u_rho_x 0.094 0.086 0.100 Latent field
15 u_rho_x 0.058 0.019 0.026 Latent field
16 u_rho_x 0.045 0.031 0.028 Latent field
17 u_rho_x 0.020 0.031 0.049 Latent field
18 u_rho_x 0.063 0.031 0.029 Latent field
19 u_rho_x 0.097 0.044 0.051 Latent field
20 u_rho_x 0.038 0.023 0.018 Latent field
21 u_rho_x 0.045 0.096 0.075 Latent field
22 u_rho_x 0.164 0.131 0.126 Latent field
23 u_rho_x 0.057 0.038 0.044 Latent field
24 u_rho_x 0.105 0.036 0.059 Latent field
25 u_rho_x 0.045 0.062 0.075 Latent field
26 u_rho_x 0.065 0.033 0.034 Latent field
27 u_rho_x 0.099 0.099 0.061 Latent field
28 u_rho_x 0.030 0.034 0.064 Latent field
29 u_rho_x 0.072 0.087 0.106 Latent field
30 u_rho_x 0.070 0.042 0.055 Latent field
31 u_rho_x 0.067 0.048 0.063 Latent field
32 u_rho_x 0.035 0.092 0.102 Latent field
1 us_rho_x 0.134 0.170 0.141 Latent field
2 us_rho_x 0.135 0.107 0.099 Latent field
3 us_rho_x 0.106 0.084 0.101 Latent field
4 us_rho_x 0.105 0.082 0.101 Latent field
5 us_rho_x 0.079 0.067 0.076 Latent field
6 us_rho_x 0.203 0.194 0.133 Latent field
7 us_rho_x 0.057 0.061 0.043 Latent field
8 us_rho_x 0.075 0.068 0.099 Latent field
9 us_rho_x 0.070 0.055 0.062 Latent field
10 us_rho_x 0.062 0.073 0.051 Latent field
11 us_rho_x 0.088 0.085 0.089 Latent field
12 us_rho_x 0.109 0.087 0.103 Latent field
13 us_rho_x 0.073 0.054 0.078 Latent field
14 us_rho_x 0.057 0.085 0.061 Latent field
15 us_rho_x 0.057 0.070 0.047 Latent field
16 us_rho_x 0.062 0.081 0.083 Latent field
17 us_rho_x 0.040 0.058 0.091 Latent field
18 us_rho_x 0.110 0.114 0.091 Latent field
19 us_rho_x 0.062 0.127 0.064 Latent field
20 us_rho_x 0.106 0.138 0.097 Latent field
21 us_rho_x 0.096 0.085 0.095 Latent field
22 us_rho_x 0.114 0.075 0.102 Latent field
23 us_rho_x 0.092 0.072 0.088 Latent field
24 us_rho_x 0.035 0.037 0.054 Latent field
25 us_rho_x 0.062 0.098 0.087 Latent field
26 us_rho_x 0.079 0.080 0.076 Latent field
27 us_rho_x 0.052 0.031 0.029 Latent field
28 us_rho_x 0.059 0.043 0.064 Latent field
29 us_rho_x 0.072 0.108 0.095 Latent field
30 us_rho_x 0.084 0.066 0.075 Latent field
31 us_rho_x 0.051 0.058 0.054 Latent field
32 us_rho_x 0.109 0.040 0.081 Latent field
1 u_rho_xs 0.128 0.080 0.081 Latent field
2 u_rho_xs 0.087 0.111 0.089 Latent field
3 u_rho_xs 0.112 0.123 0.133 Latent field
4 u_rho_xs 0.119 0.104 0.159 Latent field
5 u_rho_xs 0.193 0.163 0.152 Latent field
6 u_rho_xs 0.268 0.227 0.227 Latent field
7 u_rho_xs 0.190 0.200 0.176 Latent field
8 u_rho_xs 0.107 0.107 0.114 Latent field
9 u_rho_xs 0.123 0.139 0.135 Latent field
10 u_rho_xs 0.084 0.078 0.090 Latent field
11 u_rho_xs 0.105 0.098 0.103 Latent field
12 u_rho_xs 0.246 0.204 0.225 Latent field
13 u_rho_xs 0.104 0.135 0.129 Latent field
14 u_rho_xs 0.138 0.136 0.146 Latent field
15 u_rho_xs 0.159 0.151 0.149 Latent field
16 u_rho_xs 0.103 0.104 0.085 Latent field
17 u_rho_xs 0.102 0.095 0.104 Latent field
18 u_rho_xs 0.187 0.145 0.192 Latent field
19 u_rho_xs 0.172 0.150 0.152 Latent field
20 u_rho_xs 0.130 0.101 0.116 Latent field
21 u_rho_xs 0.110 0.127 0.134 Latent field
22 u_rho_xs 0.262 0.234 0.210 Latent field
23 u_rho_xs 0.165 0.167 0.156 Latent field
24 u_rho_xs 0.211 0.185 0.190 Latent field
25 u_rho_xs 0.115 0.100 0.107 Latent field
26 u_rho_xs 0.146 0.146 0.140 Latent field
27 u_rho_xs 0.205 0.177 0.181 Latent field
28 u_rho_xs 0.139 0.118 0.139 Latent field
29 u_rho_xs 0.203 0.173 0.165 Latent field
30 u_rho_xs 0.144 0.120 0.148 Latent field
31 u_rho_xs 0.155 0.135 0.137 Latent field
32 u_rho_xs 0.135 0.124 0.116 Latent field
1 us_rho_xs 0.036 0.036 0.037 Latent field
2 us_rho_xs 0.023 0.024 0.049 Latent field
3 us_rho_xs 0.027 0.031 0.042 Latent field
4 us_rho_xs 0.016 0.035 0.053 Latent field
5 us_rho_xs 0.044 0.045 0.048 Latent field
6 us_rho_xs 0.034 0.056 0.038 Latent field
7 us_rho_xs 0.077 0.060 0.033 Latent field
8 us_rho_xs 0.035 0.049 0.040 Latent field
9 us_rho_xs 0.028 0.039 0.028 Latent field
10 us_rho_xs 0.051 0.024 0.033 Latent field
11 us_rho_xs 0.026 0.037 0.033 Latent field
12 us_rho_xs 0.052 0.050 0.031 Latent field
13 us_rho_xs 0.040 0.039 0.022 Latent field
14 us_rho_xs 0.022 0.028 0.037 Latent field
15 us_rho_xs 0.032 0.056 0.035 Latent field
16 us_rho_xs 0.021 0.050 0.030 Latent field
17 us_rho_xs 0.027 0.025 0.019 Latent field
18 us_rho_xs 0.030 0.037 0.020 Latent field
19 us_rho_xs 0.041 0.044 0.019 Latent field
20 us_rho_xs 0.030 0.044 0.022 Latent field
21 us_rho_xs 0.041 0.027 0.025 Latent field
22 us_rho_xs 0.040 0.039 0.058 Latent field
23 us_rho_xs 0.045 0.042 0.063 Latent field
24 us_rho_xs 0.038 0.038 0.046 Latent field
25 us_rho_xs 0.076 0.049 0.068 Latent field
26 us_rho_xs 0.043 0.040 0.047 Latent field
27 us_rho_xs 0.047 0.022 0.030 Latent field
28 us_rho_xs 0.041 0.036 0.038 Latent field
29 us_rho_xs 0.031 0.051 0.037 Latent field
30 us_rho_xs 0.053 0.032 0.032 Latent field
31 us_rho_xs 0.042 0.045 0.048 Latent field
32 us_rho_xs 0.031 0.031 0.035 Latent field
1 u_rho_a 0.075 0.077 0.078 Latent field
2 u_rho_a 0.074 0.084 0.074 Latent field
3 u_rho_a 0.069 0.084 0.079 Latent field
4 u_rho_a 0.081 0.082 0.074 Latent field
5 u_rho_a 0.069 0.084 0.073 Latent field
6 u_rho_a 0.086 0.082 0.075 Latent field
7 u_rho_a 0.070 0.090 0.077 Latent field
8 u_rho_a 0.073 0.080 0.075 Latent field
9 u_rho_a 0.081 0.078 0.072 Latent field
10 u_rho_a 0.080 0.087 0.072 Latent field
1 u_rho_as 0.072 0.041 0.075 Latent field
2 u_rho_as 0.072 0.063 0.077 Latent field
3 u_rho_as 0.072 0.059 0.074 Latent field
4 u_rho_as 0.078 0.064 0.072 Latent field
5 u_rho_as 0.065 0.061 0.074 Latent field
6 u_rho_as 0.069 0.060 0.073 Latent field
7 u_rho_as 0.071 0.056 0.071 Latent field
8 u_rho_as 0.060 0.061 0.073 Latent field
9 u_rho_as 0.072 0.080 0.074 Latent field
10 u_rho_as 0.078 0.085 0.076 Latent field
1 logit_phi_alpha_x 0.509 0.509 0.101 Hyper
1 log_sigma_alpha_x 0.549 0.549 0.148 Hyper
1 logit_phi_alpha_xs 0.568 0.568 0.128 Hyper
1 log_sigma_alpha_xs 0.513 0.513 0.133 Hyper
1 logit_phi_alpha_a 0.553 0.553 0.099 Hyper
1 log_sigma_alpha_a 0.580 0.580 0.122 Hyper
1 logit_phi_alpha_as 0.537 0.537 0.058 Hyper
1 log_sigma_alpha_as 0.533 0.533 0.092 Hyper
1 log_sigma_alpha_xa 0.656 0.656 0.171 Hyper
1 u_alpha_x 0.122 0.120 0.145 Latent field
2 u_alpha_x 0.126 0.111 0.140 Latent field
3 u_alpha_x 0.114 0.095 0.129 Latent field
4 u_alpha_x 0.142 0.137 0.158 Latent field
5 u_alpha_x 0.139 0.109 0.135 Latent field
6 u_alpha_x 0.104 0.083 0.094 Latent field
7 u_alpha_x 0.107 0.059 0.059 Latent field
8 u_alpha_x 0.177 0.135 0.141 Latent field
9 u_alpha_x 0.161 0.129 0.145 Latent field
10 u_alpha_x 0.185 0.132 0.136 Latent field
11 u_alpha_x 0.154 0.113 0.150 Latent field
12 u_alpha_x 0.116 0.056 0.074 Latent field
13 u_alpha_x 0.158 0.093 0.136 Latent field
14 u_alpha_x 0.087 0.047 0.072 Latent field
15 u_alpha_x 0.099 0.073 0.085 Latent field
16 u_alpha_x 0.059 0.058 0.068 Latent field
17 u_alpha_x 0.090 0.092 0.119 Latent field
18 u_alpha_x 0.123 0.140 0.119 Latent field
19 u_alpha_x 0.112 0.117 0.118 Latent field
20 u_alpha_x 0.118 0.109 0.112 Latent field
21 u_alpha_x 0.171 0.148 0.150 Latent field
22 u_alpha_x 0.119 0.116 0.150 Latent field
23 u_alpha_x 0.134 0.119 0.142 Latent field
24 u_alpha_x 0.133 0.081 0.105 Latent field
25 u_alpha_x 0.099 0.089 0.137 Latent field
26 u_alpha_x 0.117 0.071 0.101 Latent field
27 u_alpha_x 0.131 0.132 0.145 Latent field
28 u_alpha_x 0.091 0.078 0.087 Latent field
29 u_alpha_x 0.098 0.097 0.107 Latent field
30 u_alpha_x 0.109 0.087 0.119 Latent field
31 u_alpha_x 0.086 0.052 0.081 Latent field
32 u_alpha_x 0.069 0.061 0.070 Latent field
1 us_alpha_x 0.069 0.076 0.074 Latent field
2 us_alpha_x 0.070 0.062 0.053 Latent field
3 us_alpha_x 0.085 0.066 0.074 Latent field
4 us_alpha_x 0.066 0.068 0.098 Latent field
5 us_alpha_x 0.107 0.067 0.100 Latent field
6 us_alpha_x 0.094 0.063 0.078 Latent field
7 us_alpha_x 0.025 0.026 0.031 Latent field
8 us_alpha_x 0.144 0.070 0.119 Latent field
9 us_alpha_x 0.135 0.075 0.100 Latent field
10 us_alpha_x 0.131 0.067 0.098 Latent field
11 us_alpha_x 0.103 0.057 0.099 Latent field
12 us_alpha_x 0.098 0.048 0.081 Latent field
13 us_alpha_x 0.103 0.061 0.096 Latent field
14 us_alpha_x 0.064 0.051 0.024 Latent field
15 us_alpha_x 0.087 0.050 0.074 Latent field
16 us_alpha_x 0.043 0.035 0.035 Latent field
17 us_alpha_x 0.092 0.065 0.081 Latent field
18 us_alpha_x 0.084 0.076 0.081 Latent field
19 us_alpha_x 0.084 0.073 0.083 Latent field
20 us_alpha_x 0.083 0.074 0.103 Latent field
21 us_alpha_x 0.131 0.083 0.118 Latent field
22 us_alpha_x 0.184 0.106 0.137 Latent field
23 us_alpha_x 0.084 0.077 0.084 Latent field
24 us_alpha_x 0.103 0.064 0.072 Latent field
25 us_alpha_x 0.102 0.061 0.099 Latent field
26 us_alpha_x 0.074 0.061 0.075 Latent field
27 us_alpha_x 0.077 0.057 0.083 Latent field
28 us_alpha_x 0.066 0.047 0.060 Latent field
29 us_alpha_x 0.106 0.086 0.087 Latent field
30 us_alpha_x 0.072 0.040 0.074 Latent field
31 us_alpha_x 0.051 0.035 0.040 Latent field
32 us_alpha_x 0.041 0.030 0.029 Latent field
1 u_alpha_xs 0.034 0.069 0.050 Latent field
2 u_alpha_xs 0.069 0.066 0.094 Latent field
3 u_alpha_xs 0.034 0.073 0.042 Latent field
4 u_alpha_xs 0.132 0.138 0.149 Latent field
5 u_alpha_xs 0.128 0.050 0.090 Latent field
6 u_alpha_xs 0.033 0.043 0.056 Latent field
7 u_alpha_xs 0.051 0.057 0.063 Latent field
8 u_alpha_xs 0.068 0.067 0.056 Latent field
9 u_alpha_xs 0.085 0.059 0.058 Latent field
10 u_alpha_xs 0.096 0.084 0.126 Latent field
11 u_alpha_xs 0.129 0.080 0.099 Latent field
12 u_alpha_xs 0.107 0.094 0.104 Latent field
13 u_alpha_xs 0.104 0.043 0.109 Latent field
14 u_alpha_xs 0.049 0.076 0.050 Latent field
15 u_alpha_xs 0.106 0.103 0.111 Latent field
16 u_alpha_xs 0.123 0.103 0.074 Latent field
17 u_alpha_xs 0.085 0.063 0.083 Latent field
18 u_alpha_xs 0.148 0.161 0.061 Latent field
19 u_alpha_xs 0.059 0.053 0.074 Latent field
20 u_alpha_xs 0.047 0.051 0.099 Latent field
21 u_alpha_xs 0.086 0.079 0.046 Latent field
22 u_alpha_xs 0.166 0.120 0.107 Latent field
23 u_alpha_xs 0.089 0.070 0.104 Latent field
24 u_alpha_xs 0.133 0.100 0.105 Latent field
25 u_alpha_xs 0.155 0.101 0.159 Latent field
26 u_alpha_xs 0.064 0.052 0.059 Latent field
27 u_alpha_xs 0.119 0.126 0.120 Latent field
28 u_alpha_xs 0.083 0.036 0.073 Latent field
29 u_alpha_xs 0.173 0.121 0.086 Latent field
30 u_alpha_xs 0.095 0.078 0.090 Latent field
31 u_alpha_xs 0.117 0.054 0.080 Latent field
32 u_alpha_xs 0.119 0.065 0.098 Latent field
1 us_alpha_xs 0.042 0.038 0.048 Latent field
2 us_alpha_xs 0.032 0.068 0.054 Latent field
3 us_alpha_xs 0.049 0.034 0.056 Latent field
4 us_alpha_xs 0.058 0.070 0.091 Latent field
5 us_alpha_xs 0.070 0.070 0.068 Latent field
6 us_alpha_xs 0.061 0.071 0.070 Latent field
7 us_alpha_xs 0.034 0.032 0.037 Latent field
8 us_alpha_xs 0.093 0.057 0.087 Latent field
9 us_alpha_xs 0.100 0.065 0.098 Latent field
10 us_alpha_xs 0.114 0.070 0.083 Latent field
11 us_alpha_xs 0.101 0.091 0.102 Latent field
12 us_alpha_xs 0.098 0.089 0.109 Latent field
13 us_alpha_xs 0.116 0.053 0.100 Latent field
14 us_alpha_xs 0.069 0.033 0.076 Latent field
15 us_alpha_xs 0.110 0.091 0.131 Latent field
16 us_alpha_xs 0.142 0.061 0.127 Latent field
17 us_alpha_xs 0.082 0.024 0.066 Latent field
18 us_alpha_xs 0.111 0.027 0.083 Latent field
19 us_alpha_xs 0.037 0.056 0.032 Latent field
20 us_alpha_xs 0.034 0.045 0.037 Latent field
21 us_alpha_xs 0.070 0.026 0.057 Latent field
22 us_alpha_xs 0.096 0.090 0.087 Latent field
23 us_alpha_xs 0.125 0.064 0.106 Latent field
24 us_alpha_xs 0.150 0.082 0.120 Latent field
25 us_alpha_xs 0.084 0.070 0.091 Latent field
26 us_alpha_xs 0.155 0.084 0.124 Latent field
27 us_alpha_xs 0.141 0.079 0.132 Latent field
28 us_alpha_xs 0.115 0.073 0.081 Latent field
29 us_alpha_xs 0.137 0.079 0.126 Latent field
30 us_alpha_xs 0.166 0.085 0.159 Latent field
31 us_alpha_xs 0.130 0.072 0.125 Latent field
32 us_alpha_xs 0.140 0.070 0.126 Latent field
1 u_alpha_a 0.055 0.033 0.062 Latent field
2 u_alpha_a 0.069 0.036 0.064 Latent field
3 u_alpha_a 0.066 0.031 0.061 Latent field
4 u_alpha_a 0.052 0.035 0.068 Latent field
5 u_alpha_a 0.073 0.037 0.069 Latent field
6 u_alpha_a 0.060 0.025 0.065 Latent field
7 u_alpha_a 0.069 0.023 0.066 Latent field
8 u_alpha_a 0.068 0.030 0.066 Latent field
9 u_alpha_a 0.065 0.042 0.066 Latent field
10 u_alpha_a 0.064 0.036 0.065 Latent field
11 u_alpha_a 0.075 0.034 0.065 Latent field
12 u_alpha_a 0.067 0.044 0.068 Latent field
13 u_alpha_a 0.066 0.042 0.063 Latent field
1 u_alpha_as 0.067 0.094 0.077 Latent field
2 u_alpha_as 0.100 0.103 0.098 Latent field
3 u_alpha_as 0.120 0.104 0.113 Latent field
4 u_alpha_as 0.129 0.104 0.119 Latent field
5 u_alpha_as 0.155 0.103 0.103 Latent field
6 u_alpha_as 0.094 0.074 0.091 Latent field
7 u_alpha_as 0.077 0.067 0.059 Latent field
8 u_alpha_as 0.060 0.057 0.054 Latent field
9 u_alpha_as 0.059 0.040 0.062 Latent field
10 u_alpha_as 0.064 0.044 0.061 Latent field
1 u_alpha_xa 0.033 0.056 0.052 Latent field
2 u_alpha_xa 0.067 0.039 0.022 Latent field
3 u_alpha_xa 0.058 0.036 0.025 Latent field
4 u_alpha_xa 0.080 0.063 0.071 Latent field
5 u_alpha_xa 0.109 0.074 0.094 Latent field
6 u_alpha_xa 0.109 0.122 0.085 Latent field
7 u_alpha_xa 0.081 0.088 0.061 Latent field
8 u_alpha_xa 0.061 0.048 0.060 Latent field
9 u_alpha_xa 0.069 0.037 0.044 Latent field
10 u_alpha_xa 0.103 0.092 0.109 Latent field
11 u_alpha_xa 0.094 0.086 0.139 Latent field
12 u_alpha_xa 0.116 0.068 0.100 Latent field
13 u_alpha_xa 0.073 0.017 0.065 Latent field
14 u_alpha_xa 0.084 0.094 0.075 Latent field
15 u_alpha_xa 0.054 0.043 0.050 Latent field
16 u_alpha_xa 0.038 0.052 0.042 Latent field
17 u_alpha_xa 0.046 0.032 0.024 Latent field
18 u_alpha_xa 0.038 0.037 0.025 Latent field
19 u_alpha_xa 0.098 0.053 0.083 Latent field
20 u_alpha_xa 0.042 0.039 0.059 Latent field
21 u_alpha_xa 0.082 0.077 0.078 Latent field
22 u_alpha_xa 0.218 0.184 0.188 Latent field
23 u_alpha_xa 0.024 0.037 0.055 Latent field
24 u_alpha_xa 0.061 0.020 0.050 Latent field
25 u_alpha_xa 0.091 0.063 0.093 Latent field
26 u_alpha_xa 0.033 0.053 0.033 Latent field
27 u_alpha_xa 0.097 0.122 0.090 Latent field
28 u_alpha_xa 0.069 0.045 0.054 Latent field
29 u_alpha_xa 0.035 0.050 0.022 Latent field
30 u_alpha_xa 0.058 0.062 0.041 Latent field
31 u_alpha_xa 0.026 0.035 0.019 Latent field
32 u_alpha_xa 0.038 0.027 0.042 Latent field
1 OmegaT_raw 0.512 0.512 0.024 Hyper
1 log_betaT 0.681 0.681 0.247 Hyper
1 log_sigma_lambda_x 0.630 0.630 0.175 Hyper
1 ui_lambda_x 0.112 0.124 0.097 Latent field
2 ui_lambda_x 0.123 0.121 0.110 Latent field
3 ui_lambda_x 0.101 0.126 0.105 Latent field
4 ui_lambda_x 0.108 0.121 0.107 Latent field
5 ui_lambda_x 0.124 0.101 0.134 Latent field
6 ui_lambda_x 0.093 0.092 0.104 Latent field
7 ui_lambda_x 0.114 0.136 0.124 Latent field
8 ui_lambda_x 0.115 0.119 0.116 Latent field
9 ui_lambda_x 0.117 0.097 0.127 Latent field
10 ui_lambda_x 0.175 0.165 0.139 Latent field
11 ui_lambda_x 0.121 0.108 0.130 Latent field
12 ui_lambda_x 0.110 0.088 0.109 Latent field
13 ui_lambda_x 0.117 0.118 0.121 Latent field
14 ui_lambda_x 0.114 0.104 0.113 Latent field
15 ui_lambda_x 0.103 0.137 0.108 Latent field
16 ui_lambda_x 0.114 0.097 0.134 Latent field
17 ui_lambda_x 0.124 0.103 0.142 Latent field
18 ui_lambda_x 0.105 0.101 0.097 Latent field
19 ui_lambda_x 0.093 0.135 0.136 Latent field
20 ui_lambda_x 0.122 0.113 0.107 Latent field
21 ui_lambda_x 0.120 0.112 0.105 Latent field
22 ui_lambda_x 0.098 0.104 0.121 Latent field
23 ui_lambda_x 0.117 0.105 0.094 Latent field
24 ui_lambda_x 0.104 0.126 0.109 Latent field
25 ui_lambda_x 0.125 0.108 0.104 Latent field
26 ui_lambda_x 0.143 0.126 0.145 Latent field
27 ui_lambda_x 0.122 0.101 0.098 Latent field
28 ui_lambda_x 0.087 0.093 0.098 Latent field
29 ui_lambda_x 0.079 0.112 0.103 Latent field
30 ui_lambda_x 0.116 0.117 0.112 Latent field
31 ui_lambda_x 0.123 0.122 0.101 Latent field
32 ui_lambda_x 0.127 0.099 0.130 Latent field
1 log_sigma_ancrho_x 0.502 0.502 0.027 Hyper
1 log_sigma_ancalpha_x 0.507 0.507 0.110 Hyper
1 ui_anc_rho_x 0.046 0.077 0.063 Latent field
2 ui_anc_rho_x 0.030 0.080 0.046 Latent field
3 ui_anc_rho_x 0.059 0.056 0.081 Latent field
4 ui_anc_rho_x 0.050 0.089 0.066 Latent field
5 ui_anc_rho_x 0.059 0.059 0.023 Latent field
6 ui_anc_rho_x 0.082 0.030 0.055 Latent field
7 ui_anc_rho_x 0.029 0.051 0.082 Latent field
8 ui_anc_rho_x 0.017 0.049 0.023 Latent field
9 ui_anc_rho_x 0.047 0.040 0.039 Latent field
10 ui_anc_rho_x 0.027 0.115 0.115 Latent field
11 ui_anc_rho_x 0.025 0.054 0.053 Latent field
12 ui_anc_rho_x 0.108 0.035 0.023 Latent field
13 ui_anc_rho_x 0.025 0.024 0.034 Latent field
14 ui_anc_rho_x 0.043 0.067 0.080 Latent field
15 ui_anc_rho_x 0.045 0.084 0.070 Latent field
16 ui_anc_rho_x 0.062 0.027 0.025 Latent field
17 ui_anc_rho_x 0.078 0.041 0.038 Latent field
18 ui_anc_rho_x 0.045 0.043 0.086 Latent field
19 ui_anc_rho_x 0.024 0.054 0.037 Latent field
20 ui_anc_rho_x 0.046 0.047 0.076 Latent field
21 ui_anc_rho_x 0.046 0.075 0.050 Latent field
22 ui_anc_rho_x 0.106 0.082 0.077 Latent field
23 ui_anc_rho_x 0.033 0.066 0.044 Latent field
24 ui_anc_rho_x 0.055 0.097 0.084 Latent field
25 ui_anc_rho_x 0.037 0.085 0.093 Latent field
26 ui_anc_rho_x 0.030 0.066 0.059 Latent field
27 ui_anc_rho_x 0.077 0.112 0.019 Latent field
28 ui_anc_rho_x 0.026 0.032 0.016 Latent field
29 ui_anc_rho_x 0.116 0.048 0.030 Latent field
30 ui_anc_rho_x 0.035 0.028 0.030 Latent field
31 ui_anc_rho_x 0.027 0.036 0.042 Latent field
32 ui_anc_rho_x 0.025 0.051 0.079 Latent field
1 ui_anc_alpha_x 0.074 0.064 0.072 Latent field
2 ui_anc_alpha_x 0.083 0.062 0.039 Latent field
3 ui_anc_alpha_x 0.115 0.137 0.129 Latent field
4 ui_anc_alpha_x 0.074 0.063 0.045 Latent field
5 ui_anc_alpha_x 0.036 0.071 0.045 Latent field
6 ui_anc_alpha_x 0.038 0.054 0.044 Latent field
7 ui_anc_alpha_x 0.079 0.069 0.053 Latent field
8 ui_anc_alpha_x 0.051 0.057 0.070 Latent field
9 ui_anc_alpha_x 0.091 0.071 0.092 Latent field
10 ui_anc_alpha_x 0.049 0.039 0.039 Latent field
11 ui_anc_alpha_x 0.050 0.039 0.043 Latent field
12 ui_anc_alpha_x 0.044 0.076 0.081 Latent field
13 ui_anc_alpha_x 0.042 0.063 0.046 Latent field
14 ui_anc_alpha_x 0.077 0.086 0.077 Latent field
15 ui_anc_alpha_x 0.110 0.077 0.068 Latent field
16 ui_anc_alpha_x 0.106 0.124 0.060 Latent field
17 ui_anc_alpha_x 0.084 0.100 0.103 Latent field
18 ui_anc_alpha_x 0.158 0.175 0.046 Latent field
19 ui_anc_alpha_x 0.068 0.107 0.057 Latent field
20 ui_anc_alpha_x 0.077 0.051 0.060 Latent field
21 ui_anc_alpha_x 0.086 0.064 0.071 Latent field
22 ui_anc_alpha_x 0.069 0.047 0.032 Latent field
23 ui_anc_alpha_x 0.044 0.053 0.077 Latent field
24 ui_anc_alpha_x 0.055 0.036 0.051 Latent field
25 ui_anc_alpha_x 0.098 0.086 0.083 Latent field
26 ui_anc_alpha_x 0.069 0.041 0.094 Latent field
27 ui_anc_alpha_x 0.044 0.071 0.061 Latent field
28 ui_anc_alpha_x 0.077 0.060 0.057 Latent field
29 ui_anc_alpha_x 0.148 0.168 0.061 Latent field
30 ui_anc_alpha_x 0.058 0.046 0.048 Latent field
31 ui_anc_alpha_x 0.049 0.046 0.062 Latent field
32 ui_anc_alpha_x 0.066 0.087 0.065 Latent field
1 log_sigma_or_gamma 0.518 0.518 0.036 Hyper
1 log_or_gamma 0.028 0.103 0.139 Latent field
2 log_or_gamma 0.037 0.081 0.077 Latent field
3 log_or_gamma 0.051 0.026 0.023 Latent field
4 log_or_gamma 0.029 0.050 0.049 Latent field
5 log_or_gamma 0.030 0.097 0.081 Latent field
6 log_or_gamma 0.031 0.051 0.065 Latent field
7 log_or_gamma 0.026 0.023 0.022 Latent field
8 log_or_gamma 0.046 0.024 0.029 Latent field
9 log_or_gamma 0.039 0.086 0.062 Latent field
10 log_or_gamma 0.021 0.044 0.032 Latent field
11 log_or_gamma 0.043 0.017 0.023 Latent field
12 log_or_gamma 0.050 0.117 0.140 Latent field
13 log_or_gamma 0.029 0.062 0.055 Latent field
14 log_or_gamma 0.050 0.028 0.031 Latent field
15 log_or_gamma 0.035 0.042 0.050 Latent field
16 log_or_gamma 0.029 0.050 0.057 Latent field
17 log_or_gamma 0.022 0.060 0.063 Latent field
18 log_or_gamma 0.027 0.128 0.112 Latent field
19 log_or_gamma 0.028 0.045 0.063 Latent field
20 log_or_gamma 0.032 0.031 0.030 Latent field
21 log_or_gamma 0.050 0.042 0.060 Latent field
22 log_or_gamma 0.034 0.025 0.029 Latent field
23 log_or_gamma 0.037 0.049 0.063 Latent field
24 log_or_gamma 0.024 0.054 0.050 Latent field
25 log_or_gamma 0.021 0.069 0.093 Latent field
26 log_or_gamma 0.045 0.043 0.022 Latent field
27 log_or_gamma 0.034 0.020 0.046 Latent field
28 log_or_gamma 0.034 0.020 0.048 Latent field
29 log_or_gamma 0.041 0.042 0.057 Latent field
30 log_or_gamma 0.030 0.059 0.052 Latent field
31 log_or_gamma 0.039 0.069 0.084 Latent field
32 log_or_gamma 0.025 0.040 0.027 Latent field
 
ks_summary %>%
  filter(Type != "Hyper") %>%
  ks_plot_many("TMB", "aghq")
## Warning: Removed 116 rows containing missing values (`geom_tile()`).
## Picking joint bandwidth of 0.00587

(ks_summary_out <- ks_summary %>%
  group_by(Type) %>%
  summarise(
    TMB = mean(`KS(TMB, tmbstan)`),
    aghq = mean(`KS(aghq, tmbstan)`),
    adam = mean(`KS(adam, tmbstan)`)
  )) %>%
  gt::gt()
Type TMB aghq adam
Hyper 0.12440833 0.56767500 0.56767500
Latent field 0.08093662 0.07512077 0.08234904

Save some things for output:

saveRDS(ks_summary_out, "ks_summary_out.rds")

3.3 Investigation into large KS values

Want to create rank order lists of largest KS differences between methods:

ks_comparison <- function(ks_summary, method1, method2, n) {
  ks_method1 <- paste0("KS(", method1, ", tmbstan)")
  ks_method2 <- paste0("KS(", method2, ", tmbstan)")
  ks_summary[["KS difference"]] <- ks_summary[[ks_method1]] - ks_summary[[ks_method2]]
  ks_summary %>%
    filter(Type != "Hyper") %>%
    arrange(desc(`KS difference`)) %>%
    head(n = n)
}

Nodes where TMB beats aghq:

ks_comparison(ks_summary, method1 = "aghq", method2 = "TMB", 10)
histogram_and_ecdf(par = "ui_anc_alpha_x", i = 18)

histogram_and_ecdf(par = "u_alpha_xs", i = 18)

histogram_and_ecdf(par = "ui_anc_alpha_x", i = 29)

Nodes where aghq beats TMB:

ks_comparison(ks_summary, method1 = "TMB", method2 = "aghq", 10)  
histogram_and_ecdf(par = "us_alpha_xs", i = 30)

histogram_and_ecdf(par = "us_alpha_xs", i = 31)

histogram_and_ecdf(par = "us_alpha_xs", i = 26)

3.4 Correlation between KS values and ESS

Is there any correlation between the value of \(\text{KS}(\texttt{method}, \texttt{tmbstan})\) for a particular parameter and the ESS of that parameter from tmbstan output?

ks_summary %>%
  pivot_longer(
    col = starts_with("KS"),
    names_to = "Method",
    values_to = "KS"
  ) %>%
  mutate(
    Method = fct_recode(Method,
      "adam" = "KS(adam, tmbstan)",
      "aghq" = "KS(aghq, tmbstan)",
      "TMB" = "KS(TMB, tmbstan)"
    )
  ) %>%
  group_by(Parameter) %>%
  mutate(
    par_num = case_when(
      max(index) > 1 ~ paste0(Parameter, "[", index, "]"),
      TRUE ~ Parameter
    )
  ) %>%
  ungroup() %>%
  left_join(data.frame(ess) %>%
    tibble::rownames_to_column("par_num")
  ) %>%
  ggplot(aes(x = ess, y = KS)) +
    geom_point(alpha = 0.5) +
    geom_smooth(method = "lm", color = "#56B4E9") +
    facet_grid(~ Method) +
    theme_minimal() +
    labs(x = "ESS", y = "KS(method, tmbstan)")
## Joining, by = "par_num"
## `geom_smooth()` using formula = 'y ~ x'

4 MMD

#' To write!

5 PSIS

Suppose we have two sets of samples from the posterior. For each sample we are going to want to evaluate the log-likelihood, so that we can calculate the log-likelihood ratios. We can extract the TMB objective function for the log-likelihood as follows:

tmb$fit$obj$fn()
## [1] 3327.999
## attr(,"logarithm")
## [1] TRUE
#' To write!